Hydrodynamics of thermal granular convection 
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A hydrodynamic theory is formulated for buoyancy-driven ("thermal") granular convection, recently 
predicted in molecular dynamic simulations and observed in experiment. The limit of a dilute flow 
is considered. The problem is fully described by three scaled parameters. The convection occurs via 
a supercritical bifurcation, the inelasticity of the collisions being the control parameter. The theory 
is expected to be valid for small Knudsen numbers and nearly elastic grain collisions. 
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As we know from experience, hot fluid rises. Is the 
same statement true for granular fluid, where the role of 
temperature is played by fluctuations of the grain veloc- 
ities? There is strong recent evidence that the answer 
to this question is positive. Buoyancy-driven "thermal" 
granular convection was first observed in molecular dy- 
namic (MD) simulations of granular gas in two dimen- 
sions 0] , with no shear or time-dependence introduced by 
the system boundaries. It was found that thermal gran- 
ular convection appears via a supercritical bifurcation, 
with inelastic collision losses being the control parameter 
[jlj . Strong evidence for thermal granular convection was 
recently obtained in experiment with a highly fluidized 
3D granular flow Q (see also an earlier work Q ) . In these 
two systems p]|| the convection is driven by a negative 
vertical granular temperature gradient Q which makes 
this convection similar to the classical Rayleigh-Benard 
convection |5) and its analogs in compressible fluid |i]-[)) . 
In the Rayleigh-Benard problem a negative temperature 
gradient is imposed externally. In a granular flow driven 
from below, it develops spontaneously because of the in- 
elasticity of particle collisions (Tcfl . 

The phenomenon of "thermal" convection in granular 
fluids is fascinating, as it gives one more example of sim- 
ilarities /differences between the granular and "classical" 
fluids [|llj . Though basic properties of thermal granular 
convection were investigated in the MD simulations [Q, 
no theory exists yet. The objective of the present work 
is to formulate such a theory. We will work in the regime 
where the "standard" granular hydrodynamic equations 
in 2D p2 |, systematically derivable from kinetic equa- 
tions Jl3[, are expected to be accurate. As it has become 
clear by now this requires (in addition to the strong 
inequality K -C 1, see below, and sufficiently low den- 
sity) that particle collisions be nearly elastic: q -C 1, 
where q = (1 — r)/2 is the inelasticity coefficient and r 
is the normal restitution coefficient. The nearly elastic 
limit is motivated by the MD simulations H where con- 
vection was observed already at very small inelasticities: 
4xl(T 4 < q < 2xl0~ 2 . As in the MD simulations 0, we 
will assume a velocity-independent restitution coefficient. 



We will limit ourselves to dilute flow, n <C n c , where n 
is the number density of the grains and n c is the close- 
packing density. As thermal granular convection does 
not necessarily involve clustering, the latter assumption 
is not too restricting. 

Here is an outline of the rest of this Communication. 
We will see that the hydrodynamic problem of thermal 
granular convection is fully determined by three scaled 
parameters: the Froude number F, Knudsen number 
K <C 1 and inelasticity coefficient q <C 1, and by the 
aspect ratio of the system. The translationally symmet- 
ric static steady state of the system plays the role of the 
"simple conducting state" of the Rayleigh-Benard prob- 
lem. By employing the Lagrangian mass coordinate, we 
find this steady state analytically. Then, by solving the 
granular hydrodynamic equations in a square box by a 
Latticc-Boltzmann method, we observe a supercritical bi- 
furcation at a critical value of q = q c and steady convec- 
tion at q > q c , in qualitative agreement with MD simu- 
lations [G] and experiment We then investigate the 
dependence of the convection threshold q c on K and F. 
Our results open the way to a systematic investigation of 
thermal granular convection. 

Let N 3> 1 identical smooth hard disks with diame- 
ter d and mass m move without friction and inelastically 
collide inside a two-dimensional box with lateral dimen- 
sion L x and height L y . The aspect ratio of the system 
A = L x /L y . The gravity acceleration g is in the nega- 
tive y direction. The particles are driven by a base that 
is kept at temperature T . For simplicity, the three other 
walls are assumed elastic. The hydrodynamic descrip- 
tion deals with coarse-grained fields: the number density 
of grains n(r,i), granular temperature T(r,t) and mean 
velocity of grains v(r, t). The governing equations can be 
written, in the dilute limit, in the following scaled form: 



dn 
~dt 



nV-v = 0, 



ndv/dt = V-P-Fne y 
n dT/dt + nT V • v = 
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Here e y is the unit vector along y, d/dt is the total 
derivative, P 



-nTI 



ijfT^D is the stress ten- 
sor, D = (1/2) (Vv + (Vv) ) is the rate of deformation 
tensor, D = D — (1/2) tr (D ) I is the deviatoric part of 
D and I is the identity tensor. In the dilute limit the 
bulk viscosity can be neglected compared to the shear 
viscosity Q. In addition, we neglected the small vis- 
cous heating term in the heat balance equation ( |3|). The 
three scaled parameters entering Eqs. ®) and (H) are 
the Froude number F — mgLy/To, the Knudsen number 
K = 2k~ 1 / 2 (dL y , and the collision losses param- 

eter R = 8qK~ 2 . R will be used through the rest of this 
Communication instead of the inelasticity coefficient q. 
The Knudsen number K is of order of the ratio of the 
mean free path of the grains to the system height. For 
hydrodynamics to be valid we should demand K <C 1. 
The units of distance, time, velocity, density and tem- 
perature in the scaled equations are L y , L y /T^ 2 , T^ 2 , 
(n) and To, respectively. Finally, (n) — N/(L x L y ) is the 
average number density of the grains. 

The physical meaning of the scaled numbers F, K and 
R is clear. The Froude number F determines the rela- 
tive role of the maximum potential energy of grains in 
the gravity field and their maximum fluctuation energy 
supplied by the driving base. The Knudsen number K 
determines the efficiency of the momentum and energy 
transport in the system. In the hard-sphere fluid we are 
working with, the kinematic viscosity and thermal diffu- 
sivity are equal to each other, so the Prandtl number is 
equal to 1. The inelastic heat loss number R determines 
the relative role of the inelastic heat losses and heat con- 
duction. 

The boundary conditions are T(x, y = 0, t) = 1 with 
a zero heat flux at the rest of the boundaries. Also, we 
demand zero normal components of the velocity and slip 
(no-stress) conditions at all boundaries. The total num- 
ber of particles is conserved: 



1 

A 



dx / dyn(x, y,t) = 1 



(4) 



Therefore, in the hydrodynamic formulation, the problem 
is characterized by F, K and R and the aspect ratio A, 
instead of the full set of 8 parameters m, d, q, L x , L y , g, T 
and N. 

Translationally symmetric static steady states are de- 
scribed by the one-dimensional equations considered in 
many works (see, e.g. Ref. fill): 



and 



(n s T s y + Fn s = 



{Tl' 2 T' s )'~Rn 2 s T^ 2 = Q 



(5) 



(G) 



(primes denote y-derivatives) . In our problem these 
equations are complemented by the boundary conditions 
T s {y = 0) = 1 and T' s (y = 1) = and normalization con- 
dition J Q dy n s (y) = 1 . A static state is characterized by 
the scaled numbers F and R. Equations (||) and (^|) can 
be solved by going over to the Lagrangian mass coordi- 
nate p(y) = /g n s (y')dy' [|||. In view of Eq. (]|), the 
Lagrangian mass coordinate /i changes between and 1. 
First, we solve Eq. (||) for n s (fi): 



Po- Ffi 



(7) 



where po is the (as yet unknown) pressure at the thermal 
base /i = 0. Substituting this relation into Eq. (^), we 

1/2 

obtain a linear equation for Y(/i) = T s (/i): 



(A - p) Y" - Y' - (R/2) (A - n) Y = , 



(8) 



where A = po/F and the primes now stand for /it- 
derivatives. The general solution of Eq. (p[) i s a lin- 
ear combination of the Bessel functions Jo[yl/2(A — p)\ 
and Ko[y/R/2(X — n)]. The two arbitrary constants are 
found from the boundary conditions Y(p = 0) = 1 and 
Y'(ji = 1) = 0. Now we employ Eq. ^ for n s (p) and de- 
termine the Eulerian coordinate y — y(fi) by calculating 
y = J Q dfx' /n s (//). Demanding that y(fx = 1) = 1, we 
find A (and, therefore, po) which completes the solution. 
An example of static temperature and density profiles 
is shown in Fig. [I]. We found that, at F = 0.1 and 
R < 0.7, the temperature difference between the lower 
and upper plates in the static solution agrees very well 
with the MD simulations results [Q. The negative tem- 
perature gradient is clearly seen in Fig. 1. At sufficiently 
large R, a denser and heavier gas is located on top of the 
underdense gas. This destabilizing factor drives convec- 
tion. The stabilizing factors are granular viscosity and 
heat conduction. 

We investigated convective (in)stability of the static 
state by solving the time-dependent hydrodynamic equa- 
tions numerically. A lattice-Boltzmann scheme, 
previously used to study the classical Rayleigh-Benard 
convection Jl6| , was employed. The scheme give accu- 
rate results for moderate density variations which was 
the case in the parameter range of this study. Like in the 
MD simulations yj], we considered a square box: A = 1. 
The initial conditions were the following: a uniform (and 
equal to 1) temperature, zero velocity and density equal 
to 1 plus a small sinusoidal perturbation. We fixed F and 
K and varied R. The presence (absence) of convection 
in the box was measured by computing (after transients 
die out) the velocity circulation C = § v ■ dl along the 
edges of the box. In all cases a zero circulation is ob- 
served at sufficiently small R, and the flow approaches 
a static steady state. We checked that the density and 
temperature profiles of the steady state, obtained in the 
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lattice-Boltzmann simulations, agree within 1.5 % with 
the analytic solutions of Eqs. (||) and (^). Convection 
always develops, via a supercritical bifurcation, when R 
exceeds a critical value R C (F, K). Figure [2] shows steady 
convection that appears in this system after transients 
decay. Figure || shows the bifurcation diagram. The 
same type of bifurcation (supercritical bifurcation) was 
observed in the MD simulations |l| . 
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FIG. 1. One-dimensional static temperature (solid line) 
and density (dashed line) profiles for F = 0.1 and R = 0.5. 
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FIG. 2. Steady-state hydrodynamic velocity field for 
F = 0.1, K = 0.02 and R = 3.6. 

We determined the convection onsets and bifurcation 
diagrams for two values of the Froude number: F = 0.05 
and 0.1, varying the Knudsen number K between 0.01 
and 0.06. The results of this series of simulations, de- 
picted in Fig. |J, clearly show that the viscosity and heat 
conduction (both of which scale like K) are stabilizing 
factors. Also, it can be seen that stronger gravity pro- 
motes convection as expected. 
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FIG. 3. Velocity circulation C along the edge of the box 
vs. R measured in hydrodynamic simulations (points), and 
the curve 0.125 (R - 3.186) 1/2 (dashed line). In this example 
F = 0.1 and K = 0.05. 
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FIG. 4. The critical value R c for the convection onset 
vs. the Knudsen number K for F = 0.05 (circles) and 0.1 
(squares). 

Transient motions in the system were investigated by 
monitoring the maximum value of the hydrodynamic ve- 
locity in the system as a function of time. After initial 
transients decay, the dynamics depends on whether one 
is in the subcritical (R < R c ) or supercritical (R > R c ) 
range. We found that, in the supercritical range, the 
maximum velocity first increases exponentially in time 
(no "overstability" ) , and then approaches, in an oscilla- 
tory way, a constant value corresponding to a steady con- 
vection. We used the growth rates to extrapolate to the 
critical values R c for the convection onset. We also found 
that the frequency of the decaying oscillations around the 
steady convection vanishes at the convection onset. 

The next theoretical step should be the linear stability 
analysis of Eqs. (|l])-(|^) around the static solutions, and a 
detailed investigation of the hydrodynamic modes of the 
system. In the spirit of pattern formation theory 
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one should also study convection in a strip infinite in the 
lateral direction, by varying the lateral wave number of 
the perturbation. This analysis is presently under way. 

Observing thermal granular convection in experiment 
requires several conditions, some of which can be strin- 
gent. A detailed discussion of these conditions is be- 
yond the scope of this Communication. However, there 
is one crucial issue that has to be discussed. A stan- 
dard method of fluidization of granular materials (used, 
in particular, in experiment ||) is vibration of the bot- 
tom plate. There are two important necessary conditions 
for observing thermal granular convection in vibroflu- 
idized granular beds. Firstly, the frequency of vibration 
of the bottom plate should be much higher than any rel- 
evant macroscopic frequency of the granulate (like the 
frequency of the bed oscillations or inverse sound travel 
time). Secondly, the vibration amplitude should be less 
than the mean free path of the granulate near the bottom 
wall. These conditions guarantee that there is no direct 
coupling between the bottom plate vibration and collec- 
tive granular motion. Additional conditions are those of 
convection instability. Our theory gives such a condition: 
R > R c (K,F). However, we obtained this condition (a) 
in the dilute limit n <C n c , and (b) for a "thermal", 
rather than vibrating, bottom plate. Limitation (a) can 
be severe unless the experiment is done in a 2D geometry 
(spherical particles rolling on a slightly inclined smooth 
surface and driven by a vibrating wall [ fL8| ) . To what ex- 
tent is limitation (b) severe? A full quantitative answer 
to this question requires solving a similar hydrodynamic 
problem, but with a different boundary condition [ fj^]l9| 
that mimics the vibrating wall more directly. Based on 
an analogy with other driven granular systems pC|, we 
expect that the results obtained for the two boundary 
conditions will not differ too much from each other, at 
least qualitatively. 

Finally, it is possible that the concept of "thermal" 
convection can be applicable to some granular systems 
driven by shear. A recent example is the longitudinal 
vortices observed in experiment on rapid granular flow 
in a chute Mj. 

To summarize, granular hydrodynamics provide a 
proper language for the problem of thermal granular con- 
vection and open the way to a systematic investigation 
of this and related phenomena. 
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